Physical tests for random numbers in simulations.

We propose three physical tests to measure correlations in random numbers used in Monte Carlo simulations. The rst test uses autocorrelation times of certain physical quantities when the Ising model is simulated with the Wol algorithm. The second test is based on random walks, and the third on blocks of n successive numbers. We apply the tests to show that recent errors in high precision simulations using generalized feedback shift register algorithms are due to short range correlations in random number sequences. We also determine the length of these correlations. PACS numbers: 02.70.Lq, 75.40Mg, 05.40+j, 05.50.+q. The Monte Carlo (MC) simulation method is a standard technique in physical sciences [1]. The key ingredient in its successful application lies in the quality of random numbers used, which are usually produced by deterministic pseudorandom number (PRN) generator algorithms. Several tests for PRN generators have been suggested [2, 3, 4] and conducted [3, 4, 5, 6] to test PRN generators, but none of them can prove that a given generator is reliable in all applications. Sometimes the inevitable correlations in their output can lead to erroneous results, as recently demonstrated in high{ precision MC simulations for some commonly used generators combined with special simulation algorithms [7, 8, 9, 10, 11]. In this letter, our purpose is to introduce physical tests which allow more precise characterizations of correlations which may cause problems in simulations. These tests are then used to demonstrate that for some generators, there exist nontrivial local correlations which are present in rather short subsequences of random numbers, and cannot always be detected by conventional test methods [2, 3, 6]. We start by generalizing the Ising model simulations of Refs. [7, 10, 11] to show that the local correlations lead to deviations in the cluster formation process of the Wol algorithm [12] which explains the errors observed in Refs. [7, 10, 11]. As a test of random numbers, we show that integrated autocorrelation times [13] of certain physical quantities are particularly sensitive measures of local correlations. We then introduce random walk and n{block tests, and show how these tests can be 1 used to precisely measure the length of correlations in PRN sequences. In particular, for the generalized feedback shift register (GFSR) generators we show that the length of correlations is very close to the longer lag parameter of the generator. The PRN generators tested in this work include GFSR algorithms [14], which are of the form xn = xn p xn q, where is the bitwise exclusive OR operator. They are denoted by Rp; recommended values for p and q (p > q) can be found e.g. in Refs. [15]. Other generators include a linear congruential generator xn = (16807 xn 1)mod (231 1) [16] known as GGL (CONG in Ref.[7]), RAN3 [17], which is a lagged Fibonacci generator xn = (xn 55 xn 24)mod 231, and a combination generator RANMAR [18]. The GFSR generators were initialized with 32{bit integers produced by GGL. Other initialization methods including the one in Ref. [19] were also checked but the results were una ected. In Refs. [7, 10, 11] problems with GFSR generators arose when the two{ dimensional Ising model was simulated on a square lattice with the Wol algorithm [12]. We carried out analogous simulations for an Ising model of linear size L = 16 at Kc = 12 ln(1 + p2), for a variety of generators listed in Table 1. Our implementation of the single cluster search algorithm follows Ref. [20]. We calculated the energy E, the magnetic susceptibility ̂, and the (normalized) size of the ipped clusters c. We then calculated the corresponding integrated autocorrelation times E; ̂, and c, by rst calculating the autocorrelation functions C(t) = hA(t0)A(t0 + t)i hA(t0)i2 hA(t0)2i hA(t0)i2 ; (1) and then following the procedure in Ref. [13]. 2 A summary of the results in Table 1 shows that based on this test, the generators can be classi ed into two categories. For the energy hEi, for example, deviations from the exact result of hEi = 1:45312 [21] for R31, R250, R521, and RAN3 are much larger than 3 where is the standard deviation [13]. In particular, hci reveals that in the erroneous cases the average ipped cluster size is biased. Most striking, however, is the behavior of the integrated autocorrelation times. For generators, which show no signi cant deviations in hEi, ĥi, or hci, results for the 's agree well with each other. However, for R31 and R250 the integrated autocorrelation times show errors of about 8%. We thus propose these quantities as particularly sensitive measures of correlations in PRN sequences. To compare with Refs. [7, 10] we also used the autocorrelation time test to study the decimation of the output of R250, i.e. we took every kth number of the PRN sequence. For k = f3; 5; 6; 7; 9; 10; 11; 12; 24; 48g the correlations vanish in agreement with k = 5 in Ref. [7] and k = 3; 5 in Ref. [10]. On the other hand, for k = 2m with m = f0; 1; 2; 3; 4; 5; 6; 7; 8g the sequences fail. These ndings agree with the theoretical result of Golomb [22], who showed that the decimation of a maximum{length GFSR sequence by powers of two results in statistically equivalent sequences. The errors in the average cluster sizes for some of the GFSR generators suggests that there must be correlations present within the O(L2) successive PRN's which are used in the single cluster formation of the Wol algorithm. To quantify the range of correlations we propose the following random walk and n{block tests. In the random walk test [23], we consider random walks on a plane which 3 is divided into four equal blocks, each of which has an equal probability to contain the random walker after a walk of length n. The test is performed N times, and the number of occurrences in each of the four blocks is compared to the expected value of N=4, using a 2 test with three degrees of freedom. The generator fails if the 2 value exceeds 7.815 in at least two out of three independent runs. This should occur with a probability of only about 3/400. Results for a group of generators with this test are presented in Table 2 with n = 1000. They are in agreement with the autocorrelation time test. No correlations for either GGL or RANMAR were observed. R250 and R521 pass the test with k = 3, but fail with values k = f1; 2; 26g, whereas R1279 passes with all k tested. The failure of RAN3 with k = 1 is consistent with previous test results [4, 6]. It is notable, that all the failures in this test were very clear, since even the smallest 2 values exceeded 40. However, RAN3 passed the test when every second or third number was used. The main di erence between the failing generators R250 and R521 (with k = 1) and the successful ones R1279 and R4423 lies in the lag parameter p which is less than n for the former and larger than n for the latter. Thus, it is plausible that the range of correlations depends on p and q. We studied this for p with the random walk test by locating the approximate value nc, above which the generators fail. The test was performed for R31, R250, R521 and R1279 with N = 106 samples. The results for nc were 32 1, 280 5, 590 5, and 1515 5, respectively, where the error estimate is the largest distance between samples close to nc. An example of the 2 values is given in Fig. 1. These results demonstrate that the correlations are local, in the sense that they have a range very close to p. 4 In order to further quantify the nature of correlations we present the n{block test [24]. In it we take a sequence fx1; x2; : : : ; xng of uniformly distributed random numbers 0 xi < 1, whose average x is calculated. If x 1=2, we choose yi = 1; otherwise yi = 0. This is repeated N times. We then perform a 2 test on variables yi with one degree of freedom. Each test is repeated three times, and the generator fails the test with xed n if at least two out of three 2 values exceed 3.841, which should occur with a probability of about 3=400. In cases of GGL, RANMAR and RAN3, we observed no correlations up to n = 104 for N = 106. For R31, R250 and R521 we performed an iterative study by varying n. When N = 106 samples were taken, the resulting correlation lengths nc were 32 1, 267 5, and 555 5, respectively. With better statistics N = 108, we observed no change for R31, whereas the estimate for R521 reduced to 525 1, and that of R250 to 251 1. The latter value was con rmed with N = 109 also. Typical values of 2 for R250 are shown in Fig. 2, where a sharp onset of correlations at nc is visible. The results of the random walk and n{block tests together show that for the GFSR generators, the origin of the errors in the simulations presented here and in Refs. [7, 8, 9, 10, 11] must be the appearance of local correlations in the probability distribution. Moreover, our tests show that the length of correlation lies very close to the value of the longer lag parameter p. These results are in accordance with a remark in Ref. [25], and quantify the observations of Ref. [9] where it was estimated that \triple" correlations exist about 400 numbers apart for R250. It is important to realize that in our Ising simulations it is the single cluster search [20] in the Wol algorithm, 5 where O(L2) successive PRN's are used in cluster formation that makes the system especially sensitive to correlations; if GFSR generators with p L2 (q p=2) are used, the results improve considerably. The same can also be achieved with a judicious decimation (e.g. k = 3) of the sequences. In conclusion, we have presented three simple physical tests for detecting correlations in PRN sequences. We have demonstrated the quality of these tests by being able to unravel the range and nature of correlations in GFSR generators that have recently been shown to produce erroneous results in Monte Carlo simulations. In particular, we have shown that the origin of these errors lies in local correlations present in the cluster formation process of the Wol algorithm. It is remarkable that these correlations cannot be seen in either careful statistical [6] or bit level tests [4] but only with the present test methods, and other special simulation algorithms [8, 9]. We believe that together with statistical and bit level tests the physical tests presented here form a rather complete test bench which can be used to develop better generators for demanding applications. Finally, we would also like to note that we have preliminary results for applying the tests presented here for GFSR generators with four \taps" of the form xn = xn 9689 xn 471 xn 314 xn 157 [26], which is basically a 3{decimation (k = 3) of the generator xn = xn 9689 xn 471 [27]. For these generators correlations appear to be very weak as expected, but for smaller lags errors can again be found [28]. We would like to thank R. Zi for correspondence and providing us with his results prior to publication. K.K. thanks Paul Coddington for useful discussions. Tampere University of Technology and the University of Helsinki 6 have provided ample computing resources. This work has been partly sup-ported by the Academy of Finland.7 Internet addresses: [email protected], [email protected],[email protected][1] K. Binder, in Applications of the Monte Carlo Method in StatisticalPhysics, edited by K. Binder (Springer{Verlag, Berlin, 1987).[2] D. E. Knuth, The Art of Computer Programming, Volume 2: Seminu-merical Algorithms, 2nd ed. (Addison-Wesley, Reading, MA, 1981).[3] G. A. Marsaglia, in Computer Science and Statistics: The Interface,edited by L. Billard, (Elsevier, Amsterdam, 1985).[4] K. Kankaala, T. Ala{Nissila, and I. Vattulainen, Phys. Rev. E 48, 4211(1993).[5] P. L'Ecuyer, Comm. ACM 31, 742 (1988).[6] I. Vattulainen, K. Kankaala, J. Saarinen, and T. Ala{Nissila,[email protected] no. 9304008.[7] A. M. Ferrenberg, D. P. Landau, and Y. J. Wong, Phys. Rev. Lett. 69,3382 (1992).[8] P. Grassberger, J. Phys. A 26, 2769 (1993)[9] P. Grassberger, Phys. Lett. A 181, 43 (1993).[10] W. Selke, A. L. Talapov, and L. N. Shchur, JETP Lett. 58, 665 (1993).[11] P. D. Coddington, [email protected] no. 9309017.8 [12] U. Wol , Phys. Rev. Lett. 62, 361 (1989).[13] U. Wol , Phys. Lett. B 228, 379 (1989).[14] T. G. Lewis and W. H. Payne, J. Assoc. Comput. Mach. 20, 456(1973).[15] J. R. Heringa, H. W. J. Blote, and A. Compagner, Int. J. Mod. Phys.C 3, 561 (1992); S. Kirkpatrick and E. P. Stoll, J. Comput. Physics40, 517 (1981).[16] S. K. Park and K. W. Miller, Comm. ACM 31, 1192 (1988); P. A.Lewis, A. S. Goodman, and J. M. Miller, IBM Syst. J. 8, 136 (1969).[17] W. H. Press, B. P. Flannery, S. A. Tenkolsky, and W. T. Vetterling,Numerical Recipes, The Art of Scienti c Computing, Fortran version(Cambridge University Press, Cambridge, 1989).[18] F. James, Comput. Phys. Commun. 60, 329 (1990); G. Marsaglia andA. Zaman, Stat. & Prob. Lett. 8, 329 (1990).[19] M. Fushimi, Appl. Math. Lett. 2, 2, 135 (1989).[20] J.-S. Wang and R. H. Swendsen, Physica A 167, 565 (1990).[21] A. E. Ferdinand and M. E. Fisher, Phys. Rev. 185, 832 (1969).[22] S.W. Golomb, Shift Register Sequences (Holden{Day, Inc., San Fran-sisco, 1967).[23] In an unpublished work [27], R. Zi has proposed a somewhat di erentrandom walk test.9 [24] In Ref. [9], a somewhat unspeci ed \block" test was used to study therange of correlations for R(17; 5;+).[25] A. Compagner and A. Hoogland, J. Comput. Physics 71, 391 (1987).[26] R. Zi , Phys. Rev. Lett. 69, 2673 (1992).[27] R. Zi , University of Michigan Preprint, 1992 (unpublished).[28] I. Vattulainen and T. Ala{Nissila, unpublished (1994).10 Table captionsTable 1. Results of simulations for the Ising model at criticality with theWol algorithm. The number of samples is N = 107 and k is the decimationparameter. The errors shown correspond to [13]. The most erroneousresults are in boldface. See text for details.Table 2. Results of the random walk test with N = 106 samples. See textfor details.11 Figure captionsFigure 1. The 2 values for R31 (inner gure) and R250 in the random walktest as a function of walk length n. Three independent runs in both casesare denoted by di erent symbols. The horizontal lines denote 2 = 7:815.Figure 2. The 2 values for R250 in the n{block test. Curves with circlesand squares correspond to N = 108 and N = 109 samples, respectively. Inboth cases three independent runs have been performed. The horizontal linedenotes 2 = 3:841.12 Generator kqhEiĥihciÊcR3113 1.46774(7) 0.564(2) 0.5664(3) 1.233(4) 1.058(3) 0.507(2)R2501103 1.45509(7) 0.548(2) 0.5474(2) 1.333(4) 1.143(4) 0.589(4)3103 1.45302(7) 0.545(2) 0.5452(2) 1.446(6) 1.226(5) 0.628(5)R5211168 1.45379(7) 0.546(2) 0.5461(2) 1.384(5) 1.182(5) 0.604(4)R12791418 1.45312(7) 0.545(2) 0.5454(2) 1.426(5) 1.215(4) 0.622(3)R22811 1029 1.45311(7) 0.545(2) 0.5456(2) 1.439(5) 1.226(5) 0.627(5)R44231 2098 1.45303(7) 0.545(2) 0.5454(2) 1.441(5) 1.226(5) 0.624(4)R96891 4187 1.45313(7) 0.546(2) 0.5455(2) 1.444(5) 1.229(5) 0.625(4)R199371 9842 1.45294(7) 0.545(2) 0.5452(2) 1.434(5) 1.220(5) 0.624(4)R444971 21034 1.45292(7) 0.545(2) 0.5452(2) 1.434(5) 1.219(5) 0.622(2)RAN31.45254(7) 0.545(2) 0.5446(2) 1.447(5) 1.231(5) 0.630(3)GGL1.45309(7) 0.545(2) 0.5454(2) 1.436(5) 1.221(5) 0.622(4)RANMAR1.45303(7) 0.545(2) 0.5452(2) 1.443(5) 1.227(5) 0.624(4)Table 1:13 Generatorkq ResultR3113 FAILR2501,2,64 103 FAILR5211,2,64 168 FAILRAN31FAILR2503 103 PASSR5213 168 PASSR12791,2,3,64 418 PASSR44231 2098 PASSRAN32,3PASSGGLPASSRANMARPASSTable 2:14 200 220 240 260 280 300n01020304050226 28 30 32 3401020304050